!> \mainpage DEPP
!!
!! \par DIFFERENTIAL EVOLUTION PARALLEL PROGRAM
!!
!! \authors
!!
!!          Jonas Joacir Radtke
!!
!!                 E-mail: jonas.radtke@gmail.com
!!             Curriculum: http://lattes.cnpq.br/7202527344656915
!!                    URL: http://paginapessoal.utfpr.edu.br/jonas
!!
!!
!!          Guilherme Bertoldo
!!
!!                 E-mail: glbertoldo@gmail.com
!!             Curriculum: http://buscatextual.cnpq.br/buscatextual/visualizacv.do?id=H769069
!!
!!
!! \par Institution
!!          Federal University of Technology - Paraná - UTFPR
!!
!!
!! \date March, 2013.
!!
!! \version 1.0
!!
!! \par SVN CHECKOUT
!!          http://depp.googlecode.com/svn/trunk/
!!
!! \section Introduction
!!
!! This is the documentation of the program DEPP, that implements the Differential
!! Evolution (DE) algorithm for non-linear optimization within a rectangular domain.
!! This algorithm also contains an hybridization with the Response Surface Methodology
!! (RSM) that may accelerate the convergence. The program deals only with maximization.
!! Minimization can be taken into account by multiplying the fitness function by -1.
!! The code implementation is MPI-FORTRAN 90.

!> \brief Main program
program depp

   use input
   use tools
   use hybrid
   use postp
   use output

   implicit none

   include 'mpif.h'

   character(10) tcpuf  !< Formatted CPU time
   integer :: i         !< Dummy variable
   integer :: iaux      !< Dummy variable
   integer :: status(mpi_status_size)  !< mpi: vector with information and source tag


   ! Initializing MPI (from now on the code is parallel)
   call mpi_init(code)


   ! Initializing CPU time
   tcpu1 = MPI_Wtime()


   ! Initializing MPI variables
   comm = mpi_comm_world
   tag = 42


   ! Checking for MPI initialization errors
   if (code /= mpi_success) then
      write(*,*) " =====  Error in the MPI initialization. Stopping...  ====="
      call mpi_abort(comm, code, code)
   endif


   ! Getting the total number of processors
   call mpi_comm_size(comm, nproc, code)


   ! Getting the ID number of each processor
   call mpi_comm_rank(comm, iproc, code)


   ! Generating seeds for the random number subroutine
   call system_clock(count = clock)
   seed = ( iproc + 1 ) * clock
   call random_seed(put = seed)


   ! Getting the input data
   call get_parameters(  folderin,  folderout,  sname,  iarq, reload, ffit, tcpu0, &
      ext, kss, kh, kp, kcm, nu, np, ng, dif, crs, detol, xmin, xmax, xrsm, xname, &
      x, fit, pop)


   ! Checks hybridization necessary condition
   if ( kh == 1 .or. kh == 2 ) then
      ! Checks if the population size is enough to apply the RSM
      call checks_rsm_condition(kh, nu, np, nf)
   end if


   ! if ipŕoc == 0, master processor writes the parameter to a file, gets the
   !                initial population and sends information to the slaves.
   ! if iproc /= 0, slave processors: receive information from master.
   if (iproc == 0) then

      ! Writting parameters to the output file
      call write_parameters(folderout, sname, reload, ffit, ext, kss, kh, kp, kcm, &
         nu, np, ng, dif, crs, detol, xmin, xmax)


      ! If reload=0, the first population is generated randomly, otherwise
      ! the population and its fitness are read from the backup file
      if ( reload == 0 ) then
         call get_initial_population(g, nu, np, xmin, xmax, pop)
      else
         call load_backup(folderout, sname, nu, np, tcpu0, g, fit, pop)
      end if


      write(*,"(/, a, i4, a, /)") "Processing the", g, "th generation..."


      ! Sending information to slaves
      do i = 1, nproc-1
         call mpi_send(  g,     1,          mpi_integer, i, tag, comm, code)
         call mpi_send(pop, np*nu, mpi_double_precision, i, tag, comm, code)
         call mpi_send(fit,    np, mpi_double_precision, i, tag, comm, code)
      end do

   else

      ! Receiving informations from master
      call mpi_recv(  g,     1,          mpi_integer, 0, tag, comm, status, code)
      call mpi_recv(pop, np*nu, mpi_double_precision, 0, tag, comm, status, code)
      call mpi_recv(fit,    np, mpi_double_precision, 0, tag, comm, status, code)

   end if


   ! If reload = 0 or 2, the fitness of the pop. must be calculated or recalculated
   if ( reload == 0 .or. reload == 2 ) then

      ! Starting the calculation of fitness of the current population
      do ind = 1, np

         ! Selecting the processor
         iaux = mod(ind, nproc-1) + 1

         ! Only the slave processors calculate the fitness
         if (iproc == iaux) then

            x = pop(ind,:)

            ! Fitness calculation loop
            fitloop1: do

               ! Asking to the external program 'ffit' the fitness of individual 'x'
               if (ext == 0) then
                  call get_fitness(folderout, sname, ffit, ind, nu, x, xname, &
                     fit(ind), estatus)
               else
                  call get_fitness_with_extrapolation(folderout, sname, ffit, &
                     ind, nu, x, xname, fit(ind), estatus)
               end if

               ! Analysing the exit status of the external program
               select case (estatus)

                  case (0) ! Success

                     exit fitloop1

                  case (1:10) ! Failure

                     ! Failure in the calculation of fitness function. Saving informations.
                     call save_fitness_failure(nu, g, ind, folderout, sname, &
                        x, estatus)

                     ! Since this is the first generation, another individual
                     ! must be created until the exit status is 'success'
                     call get_random_individual(nu, xmin, xmax, x)

                  case default

                     write(*,*) "ERROR: external program returned an unknown status. Stopping..."

                     stop

               end select

            end do fitloop1

            ! Sending informations to master
            call mpi_send(     ind,  1,          mpi_integer, 0, tag, comm, code)
            call mpi_send(       x, nu, mpi_double_precision, 0, tag, comm, code)
            call mpi_send(fit(ind),  1, mpi_double_precision, 0, tag, comm, code)

         end if

      end do

      ! Master processor receives the population and its fitness
      if (iproc == 0) then

         ! For each individual of the population
         do i = 1, np

            ! Selecting the processor
            iaux = mod(i, nproc-1) + 1

            ! Recieving the informations from slaves
            call mpi_recv(       ind,  1,          mpi_integer, iaux, tag, comm, status, code)
            call mpi_recv(pop(ind,:), nu, mpi_double_precision, iaux, tag, comm, status, code)
            call mpi_recv(  fit(ind),  1, mpi_double_precision, iaux, tag, comm, status, code)

         end do

      end if

   end if

   call mpi_barrier(comm, code)


   ! If iproc == 0, master processor finds the best individual, calculates
   !                the convergence measure and sends informations to slaves.
   ! If iproc /= 0, slave processors receive informations from master.
   if (iproc == 0) then

      ibest = 1

      do i = 1, np

         if ( fit(i) > fit(ibest)) ibest = i

         write(*,"(a, i4, a, 10(1pe23.15, 2x))") &
            "The performance of the", i, "th individual is ", fit(i), pop(i,:)

         if (reload == 0 .or. reload == 2) then
            write(21,"(2(i12), 100(2x, 1pe23.15))") g, i, fit(i), pop(i,:)
            call flush(21)
         end if

      end do

      ! Calculation of the convergence measure
      call get_convergence(kcm, nu, np, xmin, xmax, pop, cm)

      ! Sending convergence measure to all processors
      do i = 1, nproc-1
         call mpi_send(cm, 1, mpi_double_precision, i, tag, comm, code)
      end do

      write(*,"(a,1pe23.15)") "Convergence measure: ", cm

      write(20,"(i12, 3(2x, 1pe23.15))") g, sum(fit)/np, maxval(fit), cm

      call flush(20)

   else

      ! Receiving the convergence measure from the master processor
      call mpi_recv(cm, 1, mpi_double_precision, 0, tag, comm, status, code)

   end if

   call mpi_barrier(comm, code)


   ! Starting the generations loop. This loop is mantained while the convergence
   ! measure is greater than a specified tolerance (detol) and the generation number
   ! is lower than the maximum one (ng).
   do while ( cm > detol .and. g < ng )

      g = g+1

      ! If iproc == 0, master processor sends the updated population and its fitness to
      !                slave processors and calculates the RSM individual
      !                (if hybridization was enabled).
      ! If iproc /= 0, slave processors receive the updated population and its fitness
      if (iproc == 0) then

         tcpu2 = MPI_Wtime()

         tcpu = tcpu0 + tcpu2 - tcpu1

         call convert_real_to_time_format(tcpu, tcpuf)

         write(*,"(/, a, a)") "Accumulated CPU time: ", tcpuf

         write(*,"(/, a, i4, a, /)") "Processing the", g, "th generation..."

         ! Sending the population and its fitness to slave processors
         do i = 1, nproc-1
            call mpi_send(pop, np*nu, mpi_double_precision, i, tag, comm, code)
            call mpi_send(fit,    np, mpi_double_precision, i, tag, comm, code)
         end do

         ! Hybridization with RSM (uses the data of the previous generation)
         if ( kh == 1 .or. kh == 2 ) then

            ! Returns the best estimate given by the Response Surface Methodology
            call get_rsm_optimum(nf, nu, np, pop, fit, xrsm, es)

            ! If RSM fails, a mean individual is created
            if ( es /= 0 ) then

               ! Generates a mean individual
               call get_mean_individual(nu, np, pop, xrsm)

            end if

            ! Verifies the constraints. If the individual xrsm is out of range
            ! another random individual is created.
            call get_constraints_check(nu, xmin, xmax, xrsm)

            ind = 0

            ! Asking to the external program 'ffit' the fitness of individual 'xrsm'
            if (ext == 0) then
               call get_fitness(folderout, sname, ffit, ind, nu, xrsm, xname, &
                  fitrsm, estatus)
            else
               call get_fitness_with_extrapolation(folderout, sname, ffit, &
                  ind, nu, xrsm, xname, fitrsm, estatus)
            end if

            ! Analysing the exit status of the external program
            select case (estatus)

               case (0) ! Success

               case (1:10) ! Failure

                  ! Failure in the calculation of fitness function. Saving informations.
                  call save_fitness_failure(nu, g, ind, folderout, sname, &
                     xrsm, estatus)

                  ! Since it was not possible to calculate the fitness of xrsm,
                  ! the worst result will be attributed to it
                  fitrsm = -huge(1.d0)

               case default

                  write(*,*) "ERROR: external program returned an unknown status. Stopping..."

                  stop

            end select

            write(*,"(a, i4, a, 10(1pe23.15, 2x))") &
               "The performance of the", ind, "th (RSM) ind. is ", fitrsm, xrsm

            write(21,"(2(i12), 100(2x, 1pe23.15))") g, ind, fitrsm, xrsm

         end if

      else

         ! Receiving the population and its fitness from master processor
         call mpi_recv(pop, np*nu, mpi_double_precision, 0, tag, comm, status, code)
         call mpi_recv(fit,    np, mpi_double_precision, 0, tag, comm, status, code)

      end if

      ! For each individual of the population, another individual (trial individual)
      ! is created and its fitness is calculated.
      do ind = 1, np

         ! Selecting the processor
         iaux = mod(ind, nproc-1) + 1

         ! Only slave processors calculate the fitness
         if (iproc == iaux) then

            ! Creating the trial individual x
            call get_trial_individual(ind, nu, np, kss, dif, crs, pop, fit, x)

            ! Verifying the constraints. If the individual x is out of range
            ! another random individual is created.
            call get_constraints_check(nu, xmin, xmax, x)


            ! Calculating the fitness function
            fitloop2: do

               ! Asking to the external program 'ffit' the fitness of individual 'x'
               if (ext == 0) then
                  call get_fitness(folderout, sname, ffit, ind, nu, x, xname, &
                     xfit, estatus)
               else
                  call get_fitness_with_extrapolation(folderout, sname, ffit, &
                     ind, nu, x, xname, xfit, estatus)
               end if

               ! Analysing the exit status of the external program
               select case (estatus)

                  case (0) ! Success

                     exit fitloop2

                  case (1) ! Failure

                     ! Failure in the calculation of fitness function. Saving informations.
                     call save_fitness_failure(nu, g, ind, folderout, sname, &
                        x, estatus)

                     x = pop(ind,:)

                     xfit = fit(ind)

                     exit fitloop2

                  case (2:10) ! Generate another individual

                     ! Failure in the calculation of fitness function. Saving informations.
                     call save_fitness_failure(nu, g, ind, folderout, sname, &
                        x, estatus)

                     call get_random_individual(nu, xmin, xmax, x)

                  case default

                     write(*,*) "ERROR: external program returned an unknown status. Stopping..."

                     stop

               end select

            end do fitloop2

            ! Sending informations to master processor
            call mpi_send( ind,  1,          mpi_integer, 0, tag, comm, code)
            call mpi_send(   x, nu, mpi_double_precision, 0, tag, comm, code)
            call mpi_send(xfit,  1, mpi_double_precision, 0, tag, comm, code)

         end if

      end do

      call mpi_barrier(comm, code)


      ! If iproc == 0, master processor receives the information from slave processors,
      !                compares each individual with its trial and selects the best one,
      !                calculates the convergence coefficient and sends it to slaves.
      if (iproc == 0) then

         ! For each individual of the population
         do i = 1, np

            ! Selecting the processor
            iaux = mod(i, nproc-1) + 1

            ! Recieving informations from slaves
            call mpi_recv( ind,  1,          mpi_integer, iaux, tag, comm, status, code)
            call mpi_recv(   x, nu, mpi_double_precision, iaux, tag, comm, status, code)
            call mpi_recv(xfit,  1, mpi_double_precision, iaux, tag, comm, status, code)

            write(*,"(a, i4, a, 10(1pe23.15, 2x))") &
               "The performance of the",  ind, "th individual is ", xfit, x

            write(21,"(2(i12), 100(2x, 1pe23.15))") g, ind, xfit, x

            call flush(21)

            ! Selecting the best individual
            if ( xfit >= fit(ind)) then

               pop(ind,:) = x

               fit(ind) = xfit

               if ( xfit >= fit(ibest)) ibest = ind

            end if

         end do

         ! Hybridization with RSM
         if ( kh == 1 .or. kh == 2 ) then

            ! Finds the worst individual
            irsm = minloc(fit, dim = 1)

            ! Replaces the worst individual by the RSM one (if possible)
            if ( fitrsm >= fit(irsm)) then

               pop(irsm,:) = xrsm

               fit(irsm) = fitrsm

               if ( fitrsm >= fit(ibest)) ibest = irsm

            end if

         end if

         ! Calculating of the convergence measure
         call get_convergence(kcm, nu, np, xmin, xmax, pop, cm)

         ! Sending convergence measure to slave processors
         do i = 1, nproc-1
            call mpi_send(cm, 1, mpi_double_precision, i, tag, comm, code)
         end do

         write(*,"(a,1pe23.15)") "Convergence measure: ", cm

         write(20,"(i12, 3(2x, 1pe23.15))") g, sum(fit)/np, maxval(fit), cm

         call flush(20)

         ! Calculating the ellapsed CPU time
         tcpu2 = MPI_Wtime()
         tcpu = tcpu0 + tcpu2 - tcpu1

         ! Saving backup data
         call save_backup(folderout, sname, nu, np, tcpu, g, fit, pop)

      end if

      ! If iproc /= 0, slave processors receive the convergence measure from the master processor
      if ( iproc > 0 ) then
         call mpi_recv(cm, 1, mpi_double_precision, 0, tag, comm, status, code)
      end if

      call mpi_barrier(comm, code)

   end do


   ! Post processing with SQP (solution refinement)
   if (kp == 1) then
      !!!
   end if


   ! Master processor: data post processing
   if (iproc == 0) then

      tcpu2 = MPI_Wtime()
      tcpu = tcpu0 + tcpu2 - tcpu1

      call write_output_files(folderout, sname, nu, np, ibest, g, tcpu, &
         cm, xmin, xmax, fit, pop)

   end if


   ! Finishing MPI
   call mpi_finalize(code)

end program depp
